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Abstract 


In the present paper, we precisely conduct a quantum calculus method 
for the numerical solutions of PDEs. A nonlinear Schrédinger equation is 


considered. Instead of the known classical discretization methods based 
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on the finite difference scheme, Adomian method, and third modified ver- 
sions, we consider a discretization scheme leading to subdomains according 
to q-calculus and provide an approximate solution due to a specific value 
of the parameter g. Error estimates show that q-calculus may produce effi- 
cient numerical solutions for PDEs. The q-discretization leads effectively to 


higher orders of convergence provided with faster algorithms. The numer- 


ical tests are applied to both propagation and interaction of soliton-type 


solutions. 


AMS subject classifications (2020): Primary 35C08, 35Q55; Secondary 65L80, 
81Q05. 


Keywords: NLS equation; Quantum calculus; Numerical solution; Error 


estimates. 


1 Introduction 


The present paper is devoted essentially to the development of a numeri- 
cal scheme to approximate the solution of a Nonlinear Schrédinger (NLS) 
equation in a quantum calculus framework. The aim crosses in fact the 
restriction to the resolution of such an equation and goes further to show 
that q-calculus may provide a good framework for numerical solutions of 
PDEs in general. It is well known that a major literature on the numerical 
solutions of PDEs is based on classical methods, such as finite difference, 
finite elements and volumes, Fourier analysis, and recently wavelets. See 
(7, 5, 9, 10, 11, 12, 16, 19, 28}. 

The NLS equation is strongly linked to the modeling of real physical 
phenomena, such as Newton’s laws and energy conservation in classical me- 
chanics, the behavior of dynamical systems, the description of a particle in a 
nonrelativistic setting in quantum mechanics, and so on.Therefore, the NLS 
equation attracted researchers from both theoretical and applied mathemat- 
ics and physics. See [16, 19, 22, 28). 

Originally, Schrédinger’s stated a linear form describing a moving particle 


according to the model equation 
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8n2m 


Be eV) (1) 


where w is known as a wave function, m is the particle mass, h is the Planck’s 
constant, E is the energy, and V is a potential (see [13, 20, 23, 25, 29]). 

Based upon the analogy between mechanics and optics, Schrédinger ap- 
plied a perturbation method to show an equivalence between his wave func- 
tion in mechanics and Heisenberg’s matrix. This gave rise next to the time 
dependent model 


h? 
tite, = — Ap+V(x)b—~P) in RX, (2) 


known as the cubic NLS equation. Next, different variants and forms have 
been developed and investigated by researchers in different fields (see [1, 4, 
17, 24, 26]). 

The present paper is devoted to the development of a numerical method 


based on q-calculus to approximate the solution of a reduced NLS equation 


in R written on the form 


iu, + Aut f(u) =0, (x,t) EX x (to, +00), 


u(a,to) = uo(x), rEQ, (3) 
oe (0,1) =0, (x,t) € OQ x (tp, +00). 


We consider a domain 2 in R, and to is a real parameter fixed as the initial 
time, u; is the first order partial derivative in time, A is the Laplace operator, 
and — is the outward normal derivative. Moreover, OQ is the boundary of 2. 
Also, i = u(x,t) and up = uo(x) are complex valued functions. In addition, 
f is a nonlinear function of u assumed to be at least continuous. 

In [15], the stationary solutions of problem (3) have been studied using 
direct methods issued from the equation on the whole space. See also [14]. In 
[6], a Lyapunov-Sylvester method has been applied to solve numerical NLS 
and Heat equations. 

The organization of the present work will be as follows. In section 2, 
the q-calculus essential tools will be reviewed. Section 3 is devoted to the 
presentation of our main method. The discrete quantum version of a cubic 


NLS equation will be developed with the necessary analysis of convergence, 
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stability, solvability, and consistency. Section 4 is the subject of numerical 


experimentation due to our theoretical part. We conclude afterward. 


2 Quantum calculus toolkit 


One of the interesting fields of extensions of real analysis is the so-called q- 


theory, which provides some discrete and/or some refinement of continuous 


analysis in subspaces such as R,, which is composed of the discrete grid +q”, 
n€ Z, q€ (0,1). Recall that for all 2 € R*, there exists a unique n € Z such 
that q"t! < |z| < q", which guarantees some density of the set Ry in R. 

This section aims to introduce some basic concepts of q-theory. We 
present some definitions, notations, and properties of q-derivatives, which 
will be useful later. Backgrounds on q-theory may be found in [3, 21] and 
the references therein. 


For 0 < q < 1, denote 


Ry = {+q", n€Z} and R} =R} |_){O}. 


We propose in this section to recall two basic functions that are applied 


almost everywhere in g-theory and its applications. See, for example, [18]. 


Definition 1. The q-derivative of a function is defined by 


fle) lar) 
Difte)=4 =ae * 77% 
f'(0) , else, 


provided that f is differentiable at 0. 


The operator Dg is the q-analogue of the classical derivative, as indeed, 


if f is differentiable, we get 


lim D,f (x) = wi). 


q—1 dx 


Many concepts of derivatives and integration rules have been extended for 
the case of q-calculus. 
The only drawback of the q-calculus is the fact that they remain applied 


and investigated especially in harmonic functional analysis for the major part 


Tran. J. Numer. Anal. Optim., Vol. 14, No. 2, 2024, pp 330-346 


Arfaoui and Ben Mabrouk 334 


of the literature. A first step ahead has been conducted by Koornwinder 
and Swarttouw when studying Jackson’s third g-Bessel function. Their work 
motivates researchers to develop different q-differential operators. Recently, 
g-calculus returns to take place in PDEs. Indeed, consider, for example, an 
elliptic equation 

Au+ f(u, 2) =0, (4) 


where f is a suitable function, generally nonlinear in u. We may search for 
a numerical g-approximation by considering a grid points in R, instead of 
finite difference/finite elements used usually. In g-theory, we already have a 


q-analog of the Laplace operator expressed as 
qu(q-tx) — (1+ q)u(x) + u(qz) 


A,u(x) = = : 


For a = q” in R7, we get from (4), 
Un+1 — ql + q)Un + QUn-1 = =U _ qa fis 


where u,, = u(qg”) and f, is some discretization of f(u, 7). We thus obtain a 
recursive equation permitting to compute u, recursively. More about appli- 
cations of q-calculus in partial differential equations may be found in [3]. A 


widely known example in q-theory is the Bessel type equation 


A,u(z) = —A? u(z), 
u(0) = 1, u’(0) =0, 


(A € C), which has as unique solution a modified g-Bessel function. In the 
present paper, we will exploit the q-calculus to develop numerical solutions 
of some PDEs. 


3 The discrete NLS equation 


In this section, we develop the details of our numerical method. For this aim, 
we fix Q = [0,1] and to = 0. Fix also a time step J, = (1 — q)q*, and for 
k € N, we denote t, the kth instant. For n € N, we denote x, = q”, and 


hn = (1—q)q” the nonuniform space step. Denote also u* = u(t, 2p) the 
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net function and U* its numerical approximation (the solution of the discrete 


problem). We discretize problem (3) as follows, 


Unt —UF | 24 qur = (l-pQjur Oh 


: ee a +f(Un)=0. (5) 


By setting for n,k € N, 
r= re, od. oe, 
the discrete problem (5) becomes 
UR*! = qohUR_, +(1— (1+ qoh)UR + ofUR, + FE, (6) 
where F* = il, f(U*). Now, denote for k € N, 
U* = (Un) nen; 
the infinite vector of the numerical solution at the time k. Denote also 
By =1—(1+Q)on. 
We get the following dynamical infinite matrix-vector system: 
Ue+1 = A,U* + FF, (7) 
where A, is the infinite tridiagonal matrix with elements 


1—o% ae 0 


qot BY of 0 
0 got Bs of 0 


Ax = O  D -gek oF of 0. ~ 
gon Bn On 


Infinite (especially tridiagonal) matrices are met in many fields and have been 
applied widely. Finite forms are met in PDEs, such as finite difference meth- 


ods, in numerical analysis. See, for instance, [2, 8]. In mathematics, infinite 
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tridiagonal matrices are initially related to the so-called Jacobi operators (see 
[27)). 

Note in problem (7) that a main difference with classical methods, such 
as the finite difference method, is the possibility to relax one assumption on 
boundary conditions. We only need such an assumption for one extremity of 
the domain 2. 


Now, observe that, for each k, we get 
[1 — a6)? =14 [dG > 165? = lool?, 
and similarly, 
Jon? = 1+ (1+) on? > (L+9)7|on? = (1 +4)? lon, 


which means that the matrix A; is a dominant-diagonal matrix, which guar- 
antees the solvability of our discrete scheme and leads to the following theo- 


rem. 
Theorem 1. The numerical problem (7) is uniquely solvable, for k > 2n+1. 


In terms of the classical numerical schemes, such as the finite difference, 
the assumption k > n replaces the assumption | = o(h?), where | and h are 
the time and space steps for the finite difference scheme. 

To investigate the stability of the numerical scheme, we propose to apply 
the Lyapunov criterion for stability, which states that a dynamical system 
L(Uz, Up_1,.-.) = 0 is stable in the Lyapunov sense if, for any bounded 
initial solution Up, the solution U,, remains bounded for all n > 0 uniformly 


on n. Here, we will precisely prove the following result. 


Lemma 1. The solution U* is bounded independently of k whenever the 


initial solution U° is bounded. 


Proof. We will proceed by recurrence on k. Assume firstly that ||U°|| < 7, 


for some 7 positive. It follows from the assumption k > 2n+ 1 that 


Therefore, using system (6), for k = 0, we obtain 
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Un = qo,Up_, + (1— (14+ qoP)Up + of Un, + FR. 


As U® is bounded, and the nonlinear function f is continuous, we deduce that 
U' is bounded. So, assume that U* is bounded independently of k. Using 


system (6), we get 
1+ 3q¢ 


l1—q 
Using the recurrence hypothesis and again the continuity of the nonlinear 


\Ual + Fal: 


function F’, we deduce that 


211+ 9) a 
i= v] 


k 
Ces : 


where C > 0 is a constant independent of k. 


The consistency of the proposed method is evaluated by means of the 
local truncation error arising from the discretization scheme. Assuming that 


the solution wu is sufficiently regular, we get the principal part as 


iL 1—q)hn 
L(u)(a, t) = > Utt + ( a Urner + o(le, + hn). (8) 


It is clearly observable that the truncation operator £(wu) goes to 0 as n,k 
go to oo. This yields that the quantum numerical scheme is consistent at a 


minimum order | in time and space. 


To finish with the convergence of the numerical method, we apply the 
Lax—Richtmyer equivalence criterion, which states that for consistent nu- 


merical approximations, stability and convergence are equivalent. We thus 


obtain the following lemma. 


Lemma 2. As the numerical scheme is consistent and stable, it is then 


convergent. 


Indeed, recall here that we have already proved in (8) that the used scheme 
is consistent. Next, Lemma 1 yields the stability of the scheme. Consequently, 


the Lax equivalence theorem guarantees the convergence. 
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4 Numerical implementations 


We propose in this experimental part to develop numerical examples to val- 
idate the theoretical results developed in the previous sections. We will use 


an Ly discrete norm to evaluate the error between the exact solutions and 


It = (Soar) 


for any vector (series) X = (X;) eventually in L?(C). Denote u* the net 


the numerical ones as 


function u(a,t*) and U* the numerical solution. We propose to compute the 
discrete error 


Er = max ||U* — w'll2 (9) 


on the grid (z,), n > 0. 
We take for the rest f(u) = |u|?u, which gives the original cubic NLS 


equation. We next take the classical soliton-type solution 


u(x,t) = i exp (i(5en —6t+ ”)) sech (Vale —ct)+ 0), 
@ 
where a, qs, ¢, 0 = a a, y and ¢ are some appropriate constants. For 
t fixed, this function decays exponentially as |z| — oo. It is a soliton-type 
disturbance, which travels with speed c and with a-governed amplitude. See 
[5, 6, 10, 11, 12, 20, 22, 28, 29]. 


4.1 Propagation of a single soliton 


In the first experimentation, we focus on a single-soliton-type particle. The 
computations are done for 0 <a2<1land0<t< 1. We fix the q parameter 
to many different values according to the closeness to 0 or to 1. So, let 
q€{G = i, i = 1,...,7}. We also fix the soliton parameters a = 0.01, 
ds = 1, c = 0.1, and the phase parameters y = ¢ = 0. Figures 1 and 2 
illustrate two cases of the numerical solution for the propagation of a single 


soliton issued from our quantum numerical scheme. 
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Figure 1: Propagation of a single soliton, for g = = 


Figure 2: Propagation of a single soliton, for g = 3. 

Tables 1 and 2 illustrate the error estimates between the numerical solu- 
tion and the exact one for different values of the quantum parameter q, and 
for different values of the maximal time index kK. The space grid is fixed to 
N number of points. Truncating in an order N for practice feasibility of the 


system (7), we denote TrUX, the truncated vector 
Ue SO eg 
and TrF*, the truncated vector 
TEs PE inky dels 


where the upper script 7 is for the transpose. We denote similarly, Tr AW 


the truncated matrix 
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1- oF ok 0 
got pr oy O 
k k k 
TrAN — O qoxy By ag 0 
5 k ek ok 
0 0 qo3 83 o3 O 
gon Bry 
The system (7) will be approximated by 
Truk = TrAN Truk +TrFk. 
Table 1: Error estimates for a single soliton for N = 20 
K=10 
qd 1 q2 q3 74 95 % q7 
Er 1.1710-° 1.281075 2.411074 2.5210-* 2.1810-4  3.1110-4 2.871074 
K=15 
qd 1 q2 93 4 95 % q7 
Er 2.0110-® 2.3210-® 3.011075 2.841075 2.7110-5 3.0110-5 2.7710-5 
K = 20 
qd 1 q2 93 74 95 6 q7 
Er 1.27107" 1.4410-" 1.2310-§ 2.1510-& 2.5310-® 2.011075 2.761075 
Table 2: Error estimates for a single soliton for N = 50 
K=10 
qd 1 q2 93 94 95 6 q7 
Er 117e-5 1.28e-5 241le-4 2.52e-4 2.18e€-4 3.1le—4 2.87e-—4 
K=15 
qd 1 q2 93 4 95 % q7 
Er 20le-—6 2.32e-6 3.0le—5 2.84e-5 2.7le—5 3.0le—5 2.77e-—5 
K = 20 
qd 1 q2 93 4 95 % q7 
Er 127e-7 144e-7 1.23e-6 2.15e€-6 2.58e-6 2.0le—5 2.76e€—5 
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It is notable from Tables 1 and 2 that the numerical quantum discrete 
scheme converges with good error estimates. This encourages to application 


such a discretization idea for solving more complicated problems. 


4.2 Interaction of two solitons 


We consider here two solitons traveling with the same speed but in opposite 
directions in order to obtain the interaction phenomenon. The computations 
are done as in the previous experimentation on the space domain 2, = [0, 1] 
and a time space also I; = [0,1]. We fix the q parameter as previously to 


q€{G= 4, i=1,...,7}. We also fix the soliton parameters as follows: 


e For the first soliton, we put a = 1, g, = 2, c=4, p=0, and ¢=15. 


e For the second soliton, we put a = 2,25, qg, = 2, c= —4, py = 0, and 
o=—7,5. 


As in the previous experimentation, Tables 3 and 4 illustrate the error esti- 
mates between the numerical solution and the exact one for different values 
of the quantum parameter q, and for different values of the maximal time 


index Kk. The space grid is fixed to a number JN of points. 


Table 3: Error estimates for two-interacted solitons for N = 20 


kK =10 

q q q2 93 g4 95 6 q7 

Er 2.27e—-5 2.45e€-5 2.5le-4 2.76e-4 3.14e-4 3.44e-4 4.020 -—4 
K=15 

q q q2 93 q4 95 6 q7 

Er 211le—6 2.18e-6 2.05e-5 26le—5 2.85e—5 3.12e-5 4.05e-—5 
K = 20 

qd q q2 93 g4 95 6 q7 


Er 1.92e-—7 2.02e—7 2.15e-6 2.44e-6 3.2le—6 3.32e—5 3.876e—5 
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Table 4: Error estimates for two-interacted solitons for N = 50 


kK =10 

qd 71 q2 93 q4 95 6 q7 

Er 125e-5 1.34e-—5 144e-4 2.2le—4 2.48e-—4 2.3le—4 2.15e-—4 
KkK=15 

qd 71 q2 93 q4 95 % q7 

Er 1.15e-6 1.27e-6 1.32e-—5 1.17e-5 1.95e€—5 2.02e-—5 2.38e—5 
K = 20 

qd 71 q2 q3 q4 95 % q7 


Er 1.16e-—7 1.3le—7 2.13e-6 2.23e-6 2.27e-—6 2.47e—-5 2.55e—5 


As in the previous case, we notice from Tables 3 and 4 that the numerical 
quantum discretization yielded a very close approximated solution to the 
exact one. This is clearly shown by the error estimates, where the maximum 
error is estimated by 107‘ over all the values of the quantum parameter q in 
two Tables 3 and 4. This finding motivates the use of a quantum numerical 


scheme for more general and/or complicated PDEs. 


Remark 1. We know that in the case of classical finite difference and finite 
element methods, the error is generally estimated in the power of the space 
step h = Az, which in turn is related to the size N of the mesh or, equiva- 
lently, the number of the discretization points. In these cases, we obtain an 
error to the order h® or equivalently, the order N~®, for some exponent a. 
In the new q-quantum case, we obtain an error to the order of g*%,, which has 
the form of a geometric sequence with a ratio q® in the interval (0,1), which 
therefore converges to 0 more rapidly than the previous sequence N~°. We 
then gain a quick and more precise scheme. These facts make it unnecessary 
to concretely develop the numerical comparisons, as we know in advance that 
the classical schemes will give higher error, lower speed of convergence, and 
higher running time. Mathematically, it is easy to see that for a,@ > 0, it 
holds that for all ¢ > 0, there exisits No € N, such that g?% N® < e, for 
all N > No. This means that, the q-quantum scheme, even at a lower order 
8 < qa is better than the classical schemes. In other words, the q-quantum 
scheme does not necessitate calibrating the discretization to get higher or- 


ders for convergence, consistency, and good stability, as it is done for the 
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classical finite difference and finite element methods with implicit method, 
Crank—Nicolson, centering/decentering method, higher regularity order finite 


elements bases, and so on. 


5 Conclusion 


In the present paper, the principal aim was to test the efficiency of the quan- 
tum calculus in the approximation of the solutions of PDEs. As a pro- 
totypical example, we applied qg-calculus to derive a numerical scheme for 
the well-known cubic NLS equation. As expected, the q-calculus yielded 
good approximations illustrated by low error estimates. The findings in the 
present paper make therefore good motivation to continue to exploit quan- 
tum calculus for the numerical (and also exact) solutions of different types 
of PDEs. Comparisons with other models, such as finite difference, finite 
volumes, and also wavelets as recent developments in mathematical analysis 
are fascinating and motivating future extensions. Compared to the classical 
finite difference scheme method, we may conclude theoretically that, in fact, 
the present quantum scheme is more efficient, as it is based on geometric 
sequences of time and space steps, which surely converge more rapidly than 
arithmetic discretizations. Therefore, we expect that involving or including 
hybrid schemes may induce the best results. Finally, an interesting question 
raised from the present work may be formulated as follows: Given an infinite 
matrix that is truncated in an order n. We know that at most in C, any 
truncation has at most n eigenvalues. What can we expect for the original 
linear operator defined by means of the infinite matrix? This gives rise to 
possible chaotic behavior as a future study of the present case of matrices, 


which are issued from parabolic, hyperbolic PDEs. 
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